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ABSTRACT 


Heat transfer rates for laminar, convective heat transfer in the 
entrance region between parallel plates are investigated. The hydro- 
dynamic solution due to Bodia 27 was uSed in the solution of the energy 
equation in finite difference form on a digital computer. The thermal 
boundary conditions include: constant heat input, constant wall tempera- 
ture, one wall constant temperature and one wall insulated, and constant 
but different wall temperatures on the upper and lower aie ; 

The approximate, integral methods of Siegel and Sparrow /B/, Vay 
produce results that are in close agreement with the solutions in this 
analysis for the constant heat input and constant wall temperature cases. 

The scope of the finite difference solution is limited to a narrow 
range of Prandtl numbers near unity, due to the small grid sizes required 
for convergence at small Prandtl numbers and to the overly low transfer 


rates indicated near the emtrance for high Prandtl numbers, which is a 


result of the "finite starting length." 
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TABLE OF SYMBOLS 


English Letter Symbols 


Specific heat 


Half the distance of separation of the parallel flat 
plates 


Hydraulic diameter D = 4d 
Unit conductance for convection heat transfer 


Unit thermal conductivity 








| -/Dt = Gee 

Dimensionless length parameter, ~ = ae 
e 
D 

Heat"flux per unitiarea; o/A 

; qD 
Dimensionless heat flux, 

Akt, 


Temperature 

Initial flu¥d temperature 

Mixed mean temperature 

Wall temperature 

Lower wall temperature (asymmetric case) 
Upper wall temperature (asymmetric case) 
Dimensionless temperature , Cate 
Dimensionless wall temperature, tw/t, 
Dimensionless mixed mean temperature, t /to 


Axial velocity component 





Initial axial velocity component 
Dimensionless axial velocity, u/U, 
Crosswise velocity component 

Crosswise velocity component at centerline 
Initial crosswise velocity component 


Dimensionless crosswise velocity component, V/U5 


Dimensionless length parameter, + 
Uo 


Normal coordinate 


Dimensionless normal coordinate, y/d 


Greek Letter Symbols 


g 


fh 
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Fluid density 
Fluid viscosity 


Kinematic viscosity, “@ 


Non-dimensional Groups 


N 
Uy 


PR 


R 
eq 


CD 


Local Nusselt number, —fbD. 


k 


Prandtl number, 44+*— 


kK 
du 
Reynolds number, based on half the gap width, = 
Reynolds number, based on hydraulic diampter, et 


X1 





I. INTRODUCTION 


Efficient design of compact heat exchangers such as are employed 
for gas turbine regenerators, among other applications, requires maximum 
utilization of the high transfer rates available in the entrance region. In 
this region the velocity and temperature profiles of the entering fluid are 
undergoing rapid transition from their uniform distribution at entrance to the 
fully established distributions encountered farther downstream. The well 
known Graetz solution, in which only the temperature development is con- 
sidered, does not yield adequate results when the velocity and temperature 
distributions are simultaneously developing. This is particularly true if the 
entrance region represents a considerable percentage of overall length and 
the two profiles develop at approximately the same rates. This situation is 
frequently encountered with gases where the Prandtl number is approximately 
one. 

The purpose of the present analysis is to present a finite difference 
solution for the heat transfer rates between parallel flat plates for several 
different boundary conditions corresponding to various constant wall 
temperatures and constant heat input configurations. The approach to the 
problem was suggested by Miller /i/, which is adapted from the numerical 


procedures applied to the entrance region of the circular tubes by Kays 1 2g: 





While an analysis such as the present one can hardly compete for 
simplicity with approximate treatments such as the Graetz analysis, or 
the several energy integral analyses commonly employed in engineering 
design, it does provide an exact solution for a limited number of cases 


which serves as a standard with which the accuracy of the various approxi- 


mate methods may be assessed. 





II, ANALYSIS 


A. Governing Equations 

The governing equations are the energy, momentum (Navier-Stokes) 
and continuity equations. In order to reduce the complexity and coupling 
of the equations, the following assumptions are made. 
The flow is assumed to be: 

1. steady, two-dimensional 

2. laminar 

3. incompressible 
In addition it is assumed that: 

1. thermal diffusivity (C< ) is constant 

2. convective heat transfer is large compared to radiation, axial é 


conduction, and viscous dissipation. 


The governing equations in reduced form may then be expressed as: 


Energy: 
A ot ot = sé i! 
etal 33 x 34? (1) 
Momentum: 
au +V ou = rene oP + au 
“a gx ay € 2x > yy? (2) 
Continuity: 


su, av 
ety 9 (3) 





Due to the assumption of incompressibility the above equations are 
no longer coupled. The temperature development remains dependent upon 
the velocity profile, although dependence of velocity on temperature has 
been removed. 

The boundary conditions associated with equations (1), (2), and 
(3) car flow between parallel flat plates with uniform velocity and tempera- 
ture at entrance are: 


1. At the entrance (x = 0) 


t=t =C 

0 1 
U =U, = Go 
Ve = ay 


2. At the walls (y =+ d) 


4. Applied thermal boundary conditions 
(1) symmetric cases (with respect to duct centerline) 


a. constant wall temperature 


Wy=C 





b. constant heat input 


{3 


LP x0 = 04 


(2) asymmetric:cases 
a. one wall constant temperature, one wall insulated 


t = C 
wil 5 


() 40 =a 


b. constant wall temperature 





| = C 


fF 
« 


two 


ll 
Q) 
“I 


where: CF Cy 
Solution of the energy equation requires that a hydrodynamic 


solution be first obtained. 


By Hydrodynamic Solution 


Several hydrodynamic eoluticnenten the entrance region are available. 
Solutions have been presented by Bodoia VV Schlichting TES, Schiller ie 
and Han ee The approximate solutions of Schiller and Han have the 
advantage of being analytic, and therefore, less cumbersome to use. They 
are, however, inherently less accurate than the numerical solutions presented 


by Bodoia and Schlichting. 





Bodoia's solution was obtained by evaluation of the Prandtl boundary 
layer equations in finite-difference form. Velocity profiles were obtained 
at given streamwise locations by utilizing matrix methods to simultaneously 
solve a column array of momentum expressions normal to the flow. 

A comparison of Schlichting's series solution and the Bodoia solution 
indicates that, in addition to the velocity gradient discontinuity in the 
Schlichting solution where the upstream and downstream solutions are joined, 
Schlichting’s representation of the pressure gradient, based only on the 
centerline velocity, results in a core velocity gradient which is too large 
and a wall velocity gradient which is too small VV 

Comparisons of the dimensionless pressure decrement versus length 
due to Bodoia, Schlichting and Han are presented in Fig. 2. Velocity dis- 
tribution comparisons between Bodoia and Schlichting, and Schlichting 
and Schiller, are shown in Figs. 3 and 4. 


The Bodoia solution was adopted for the present analysis. 


Ci. Energy Solution 


Equations (1) and (3) may be written in dimensionless form by 


introducing the following dimensionless variables: 


= t/t X=_V x= x/a_ U =u/u 
0 du, Re . 
Y =a V< v/u, 





In dimensionless form, the equations are then: 


Energy: 
1) aT + ReV ot == At, (4) 

Continuity: 
au +Re af 20 5) 


The following finite difference approximations are then introduced: 








or ~ OT = Txes Y le Y (6) 
aX AX OX 
Oo ee ol Tx, v2 7 T'x,y-1 (7) 
oY LOY, ALG 
Soil eli pg ae totliernccia Tx yea 
Sp AY" (8) 
aU = SU = Uxrit,¥ is Ux-1,Y (9) 
aX AX 2 AX 
ov. AV = Ney-Vxy-1 (10) 
oY AY AY 


Note that (6) is a forward difference and (10) is a backward 


difference, whereas the remaining equations are of the central difference 





type. These have been expressed in this manfner for simplicity and 
convenience. The larger truncation error normally associated with the 
forward or backward difference is not significant here since practical grid 
dimensions require that: 


AY >> AX 


and, due to the nature of the flow: 
Le) AV 
AX >? fay 
Introducing (6) through (10) into (4) and (5) and solving for the 


appropriate term yields: 


4X AX ReV 
Txsay = | OpRave * U2 AY Tx,¥-1 + 
1 - £8 | Tx, +1 + 
AX ¥ eas 
upRavt~ “yaar! '*Yt te 
Vey - Vejye a wty| Uray - Urey (12) 


which constitute the computing equations. 

Note that the Reynolds number (R e) , although appearing in the 
above equations, remains only for convenience and could be eliminated. 
It is no longer a parameter in the set of equations having been included 


in the dimensionless length variable (X). 





A mixed mean temperature (t =) may be expressed in the following 


}cpeut dA 
t 


_ 1S 
ms J cpeuda — 


way: 


Removing the constants from within the integral, and simplifying, a 


dimensionless mixed mean temperature may be written in finite difference 
form: 


(14) 
Ux Y 


Due to the no-slip condition imposed upon viscous flows, heat 
transfer to or from the fluid at the wall may be expressed in terms of 


conduction through a thin film of stagnant fluid. Thus: 


q =q/A= K( 34 Jyz0 (15) 
which must equal the heat flux due to convection, 
q =q/A- h(t, - t,) (16) 
Defining the local Nusselt number as; 
- «whd. 
Nu = = (17) 
Combining (14) and (15), 
(35) 
> 20 
= 18 
Uy ripe (18) 





A Nusselt number based on the hydraulic diameter (d) instead of 


half the gap width (d) in finite difference form becomes: 


4 ( te 7 gene) 
AY 


a - 
The hydraulic diameter was introduced in (13) to enable comparisons 
with previous solutions. The conversion is: 
D = 4d 
In addition, for the asymmetric constant wall temperature case, it 
is convenient to use a dimensionless parameter other than the Nusselt 


number. A dimensionless heat flux parameter may be defined by using 


(15) above: 


+ 


4 iA 
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D. Thermal Boundary Conditions In Finite Difference Form 


Constant Wall Temperature Symmetric Case) Since the temperature 


distribution is symmetrical with respect to the centerline, computations 
were performed in only the upper half of the duct. A coordinate system was 
adapted with the X axis on the centerline of the channel. The thermal boundary 


condition at the wall was applied as: 


10 





where N is a positive integer. 


Constant Heat Input Gymmatiic Case) With the axis system as above, 
the slope of the thermal profile was maintained constant at the wall by 
adding a constant to the calculated temperature one step from the wall, 


thus: 


where: M = ( AT/ AY),, AY 

Constant Wall Temperature with One Wall Insulated. For this case 
the coordinate system was redefined so the X axis coincided with the lower 
wall. The lower wall boundary condition was taken to be: 


wie =a 
O 


while zero slope at the upper wall was obtained by equating the wall 


temperature to the calculated temperature one step from the wall: 


two/t,, 7 Two 7 Tey - AY ° 


Constant Wall Temperature(Asymmetrio Gase).. With the axis system 


the same as for the insulated wall case, thermal boundary conditions were 


established as; 


where Ny and No» are positive integers and N, # N» 2 


11 





Noting that: 
tw1 - to = twy/to - 1 ages.) 


pL §G tw,/to - j Ny - 1 


Values of N, - ie No - 1 a&1/2, 1/3, 1/4 were investigated. 


BE. Computational Procedure 


The finite difference solution was obtained by using a CDC 1604 
digital computer to solve the equations at intervals of AX = 0.001 and 
AY= 0.1. Velocity distributions were obtained by interpolating plotted 
curves of Bodoia's tabulated results. 

Since the entire solution is dependent upon the first calculations 
at the entrance, starting values were obtained by reducing the grid size 
to4xX = 0.0001 and 4 Y = 0.05 for the fipst 20 streamwise stations from 
Po— 040 0.0020. 

The continuity equation (2) was evaluated “and the cross-wise 
velocity component {V) was introduced into the energy equation (11) at 
each grid point. The energy solution was then obtained at intervals of 

Y between centerline and the wall at each streamwise station. 

Solutions for each thermal boundary condition were obtained for 
Prandtl numbers of: 0.5, 0.7, 1.0, 126). 0%2anae uo 

Solutions for Prandtl numbers less than 0.5 are excluded by the grid 
size. This may be seen from inspection of the coefficient of the second 


temperature term in (11), which indicates the interdependence of PR,4X, 


12 





and AY. In order to obtain convergence, PR must be large enough to 
prevent this term from becoming negative. Thus low Prandtl numbers, 
approaching those of liquid metals, cannot be effectively handled by 
finite-difference methods owing to the large number of calculations 
required by the small grid. 

The computer program, written in Fortran 60, is included in the 


Appendix. 
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lil. RESULTS 


Tables I through VI contain computed and extrapolated values of 
the heat transfer rates as a function of distance from the entrance for each 
of the thermal boundary conditions for a Prandtl number of 0.7. Figures 
5 through 29 contain curves of the transfer rate versus downstream position 
and temperature and velocity profiles for Prandtl numbers of 0.5, 0.7, 1.0, 
io, 3.2 and 10. 

With the exception of the Asymmetric Constant Wall Temperature 
Case, which is expressed in terms cf the dimensionless heat flux parameter, 
qv: the heat transfer rates are preserited as local Nusselt number variation 
with the length parameter, L. Whereas, the computations were carried out 
using the dimensionless length parameter X, which is normally employed 
for hydrodynamic developments: L, based on the hydraulic diameter is 
customarily employed in heat transfer. The conversion is: 

X= 16L 

The use of all finite starting length" in the finite difference equations 
results in a finite initial value for the Nusselt number at X = 0 rather than 
the theoretically infinite value at the entrance. This causes a nearly flat 
slope near the entrance and an artificial inflection point in the Nusselt 
number versus length curve. The inflection point occurs progressively 


farther downstream with the increasing Prandt! number, thereby limiting the 





the useful span of the curves. A reduction of grid size from 4X = .001, 
AY=.1to AX= .0001, AY = .05 was observed to double the initial 
Nusselt number from 40 to 80, ard move the inflection point upstream. 
This increased the useful span. However, the approach obviously reaches 
a practical limit very rapidly due to the number of calculations and the 

- accuracy required of the hydrodynamic solution very close to the entrance. 

In general, for the three thermal boundary conditions in which the 
Nusselt number was used as a heat transfer parameter, the curves for 
Prandtl numbers of 0.5, 0.7, 1.0 and 1.6 behave as expected for values 
of L> 1074. For Prandtl numbers of 3.2 and 10, L= 1073 is the approxi- 

- mate lower limit of validity. 

Constant Heat Input. A comparison of local Nusselt numbers obtained 
in the present analysis with those reported by Siegel and Sparrow Jaye and 
Han /&/, contained in Fig. 5, indicate close agreement with the results of 
Siegel and Sparrow. Their solution employed a simplified energy equation 
in which the temperature distribution in the boundary layer was expressed 
as a series of polynomials in the transverse coordinate using the down- 
stream station as a parameter. The hydrodynamic solution employed 
Schiller's approximation. 

The comparison with Han's integro-numerical solution is not as close 


as might appear in Fig. 5 since this curve is for a Prandtl number of 0.8 


1S 





and should lie above the other two curves which are for a Prandtl number 
OiOi7 

The computed curves of local Nusselt number versus length contained 
in Fig. 6. for Prandtl numbers of 0.5, 0.7, 1.0, 1.6, 3.2, and 10,:may be 
_ approximated within 10 per cent, for 0.5< PR<$1.6, and L? 107° by the 


following; 


BR 
Se sels ower ee Ly) 


Ux 1+ .0178 (PR/j)-6 


Velocity and temperature profiles for constant heat input are contained 
in Figs.. 7 through 10. 

Constant Wall Temperature (Symmetric Case). Local values of Nusselt 
numbers for this boundary condition are compared in Fig. 11 with those of 
Sparrow /5/, which were computed using the Karman-Pohlausen method and 
Schiller's hydrodynamic solution. Again close correlation is indicated 
between the two results. 

The local Nusselt number curves of Fig. 12 may be approximately 


represented within 10 per cent for 0.5¢PR$1.6, and L? Tei by the 


following; 


N. = 7,60 + 0148 (°/L) 
Ux 1+ .0143 (PR/,).6 


Velocity and temperature profiles for this boundary condition are 


contained in Figs. 13 through 16. 
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' 
One Wall Constant Temperature - One Wall Insulated. The results 


of this computation are presented in Fig. 17. The local Nusselt number may 
be approximated within 10 per cent for 0.54 PR£1.6, and L >? 10 By; 
N, = 4.84 + _20155 (PR/) : 
= ee .0F2a6 /:) : 

The velocity and temperature distributions are illustrated in Figs. 
18 through 21. 

Constant Wall Temperature (Asymmetric Case). Dimensionless heat 
flux values for this boundary condition are presented in Figs. 22, 23 and 24 


for N}- 1 _ ratios of 1/2, 1/3, and 1/4. 


Note that the upper (hotter) wall curve is smoothly asymptotic to 
the value predicted for the steady state condition. The lower curve 
exhibits a much shallower initial slope, passes through an inflection 
point, and then, approaches a negative asymptotic value of the same 
magnitude and at the same rate as the upper wave. 

Velocity and temperature profiles for N} - 1 = 1/4 are presented 

Ho = 1 
in Figs. 25 through 28. 

The behavior of the local Nusselt number for this asymmetrical case 

is depicted in Fig. 29. The lower curve, representing the local Nusselt 


number variation for the lower (cooler) wall, does not depict the actual 


physical situation. 


17 





Both the denominator and the numerator of; 


4 ( ik ~ Tweay) 
N= 44 (18) 
x ae, het yee 
approach zero at different rates. The curve tends to positive infinity as 
the denominator approaches zero, returns from negative infinity, and 


becomes zero as the numerator goes to zero; and finally approaches the 


actual steady asymptotic value. 
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IV. CONCLUSIONS 


From the results of this study it may be concluded that: 


i 


The finite difference method, as employed in this analysis, 
yields results which are consistent with previous solutions 
for other geometries Vey, and which are asymptotic to the 
known fully developed values. The method is, however, 
limited to a small range of Prandtl numbers near unity by the 
small grid sizes necessary outside this range. For small 
Prandtl numbers the grid must be very small to obtain con- 
vergence, while for large Prandtl numbers practical grid 
dimensions result in significant inaccuracies near the 
entrance due to the "finite starting length”. 

The results of the approximate integral methods suggested 

by Siegel and Sparrow Voy for constant heat input, and by 
Sparrow ey, for constant wall temperature, Compare very 
favorably with the finite difference solutions reported here. 
This additional substantiation of the methods suggested by 
these authors is significant because of the practical utility of 
their approach,which requires relatively few calculations, is 
capable of handling a wide range of Prandtl numbers, and can 


predict heat transfer rates very cloSe to the entrance. 


rg 
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TABLE I 


CONSTANT HEAT INPUT 


PR = 0.7 


iy 

46.2 
30.0 
23.7 
18.0 
16.6 
IZ 
10.3 

8.87 
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TABLE II 
CONSTANT WALL TEMPERATURE 


(SYMMETRIC CASE) 


PR=0.7 

X ie Nu 
.0016 0001 34.0 
.0032 .0002 26.0 
0064 .0004 19.9 
.0128 - .0008 ioe 
.0160 .0010 G57 
.0320 .0020 10.5 
.0640 .0040 8.58 
.1280 .0080 7.78 
.1600 .0010 16s 


Asymptotic to 7.60 
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TABLE III 


ONE WALL CONSTANT TEMPERATURE - ONE WALL INSULATED 


23 


PR = 0.7 

X L neue Te 
.0016 .0001 33.0 1.000 
- .0032 .0002 24.0 1.000 
0064 0004 FIe2 1.000 
.0128 0008 13.0 1.000 
.0160 .0010 11.6 1.000 
0320 0020 w209 1.000 
.0640 .0040 7.20 .000 
. 1280 0080 ono O01 
. 1600 -0100 Soo yz 
-4000 50250 4.96 mrSo 
Asymptotic to- 4.84 .000 





TABLE IV 
CONSTANT WALL TEMPERATURE 


(ASYMMETRIC CASE) 


1] = = 0.7 
] 

x L ci 
.0032 0002 Zon 
00064 .0004 ies 
OAS: 0008 RS, 
.0160 .0010 10.7 
0320 0020 8.08 
.0460 .0040 6.18 
1280 . 0080 4.73 
. 1600 .0100 4.24 
» 3200 0200 2.24 
.6400 .0400 -.19 
.9600 .0600 =e 

Asymptotic to -2.00 
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TABLE VI 
CONSTANT WALL TEMPERATURE 


(ASYMMETRIC CASE) 


Ni - 1=1/4 PR = 0.7 
No - ] 
* * 
BS 2 qy Gy 
ROS Z .0002 Loe2 116.0 
.0064 .0004 16.8 82.9 
.0128 .0008 aA 57a 
.0160 .0010 10.5 50.4 
-0320 .0020 8.09 Shae 
.0640 .0040 6.17 27.5 
. 1280 .0080 4.66 20.3 
. 1600 .0100 4.06 18.4 
. 3200 .0200 1.06 bes 
.6400 .0400 -2.96 9.03 
£9600 .0600 -4,71 7 ae 
Asymptotic to -6.00 +6.00 
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APPENDIX A 


SAMPLE FORTRAN 60 PROGRAM 
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pomeam trite we Ob t- 
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NOADADNANNDNDADADADANA 


500 


40 


4] 


43 
44 
45 
523 
17 


700 
701 


ad 
710 


150 


PROGRAM HEATX23 


CONSTANT HEAT INPUT 
BODOITA VELOCITY DATA 


THIS PROGRAM USES A REDUCED GRID OF. X=.0001sY=.05 

FOR THE FIRST 20 STREAMWISE CALCULATIONS. 

GRAPH OUTPUTS OF NUSSELT NUMBER VERSUS LENGTH PARAMETER 
FOR’ PRANDTL NUMBERS OF O05 s00e7eleOslebs3e2xAND 10» 

ARE OBTAINED TO 50 PER CENT OF THE HYDRODYNAMIC 
DEVELOPMENT. PRINT OUTPUT INCLUDES NUSSELT NUMBER, 
SOPE OF THE TEMPERATURE PROFILE AT THE WALL » AND THE 
SLOPE OF THE TEMPERATURE PROFILE AT THE WALL»sAND THE MEAN 
‘IXED MEAN TEMPERATUREe IN ADDITIONs VELOCITY AND 
TEMPERATURE VALUES ARE PRINTED AT EACH GRID POINT AT 
SELECTED STREAMWISE LOCATIONS. 

DIMENSION T (30923) 9U( 3023) 9V( 30923) 9¥(23) 9X( 500) 


1sXN(500) sITITLE(12). 


READ 5O0O0s(ITITLE( I)» IT=1s6) 
READ SOOs(ITITLEC( I)» IT=7912) 
FORMAT (6A8) 

DO 750 KK=156 

GO T0(40541942543,5 a 45)9KK . 
PR=0.7 

GO TO 46 

PR=0.5 

GO TO 888 

GO TO 46 

PR=1.0 

GO TO 46 

PR=1]1.6 

GO TO 46 

PR=3e2 

GO TO 46 

PR=10. 

PRINT 523 

FORMAT (1H1) 

PRINT 17,5 PR 


FORMAT (/////10X917H PRANDTL NUMBER =9F 53) Jie 
L=] , 
GO TO 150 


READ 7Ols ((U(IsJ)s9 J=13922)s 12922) 
FORMAT (10F604) 


DO 710 I=2s22 


DO 711 J=1512 
UCT sJ)=U( 1913) 
CONTINUE 
DELX=.000] 
DELY=0.05 
R=1000. 

DO 110 L=2s21 
UTM=0.0 

UM=0.0 

XL=L 


X(LI=(XL-12¢0)*06 0001 


60 
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OG@eill J=3522 
XJ=J 
Y(L)=(XJ-3200)*06 ae 
Y(23)=1.0 
Ull»sJ)=Hale 
U(1523)=1.0 
Vi lsJ-1)=0.0 
Vil»J)=06 
T(1,J)2=1.0 
T(15J+1)=1-0 
T(1,J—-1)=1-0 
T(Lol)=T(bs5) 
U(L»s1)=U(Ls5) 
U(bL523)=0e20 
U(L52)=U(Ls4) 
Vi(i»3)=0-0 
VilL»2)=V(bLo4) | 
C THERMAL BOUNDARY CONDITIONS 
T(bs23)=T(bLs22)+205 
Ti{Ls2)=T(bs4) 
IF(L~1) 10951095108 
108 ViloJ)= VibsJ~1) + DELY/(2e¢*#DELX*R) * (UC L-1,J4)eU(L+19J))°: 
109 A= DELX/ (U(LoJ)*PR*¥DELY**2 ) 
V(tiL523)=0.0 
B= (DELX¥R¥V(LoJ)) / (2e*DELY*U(L»J)) 
, TILL J) R= (CARKB)ET( Lo Jt] dtl lem2e*A)¥T(L» Be J-1) 
ares ) Lied 1 36 LZ 
113 UTM=UTM4+e5*U(LodJ)*¥T( Lod) 
UM=UM+.5*U(L»J) 
GO TO 114 
112 UTM=UTM+U(L»sJ)*¥T(LoJ) 
UM=UM+U (LJ) 
114 TM=UTM/UM 
L111 CONTINUE 
Tl= ( T(bs23)-T(bs22)) 7 DELY 
T2= Ti{Ls»23) -~-— TM 
XN(L)=4-0*T1/T2 
131 PRINT 1185 X(L)oXN(L)sT2sT1l9TM a. 
118 FORMAT (// 4H X =9F60e4s5X914H NUSSELT NOe =9F8e595X98H TWHTM =o 
1F8e¢595X98H SLOPE =9F8e595X95H TM =9F8e5) 
PRINT1199 (TlbsJ)s J=1912) 
PRINT119» (TlLsJ)s J=13923) 
119 FORMAT (/9H TEMP. =912F704) 
PRINT1205-(U(LsJ)5 J=1912) 
PRINT1205 (U(bLoJ)»9 J=13923) 
120 FORMAT (/ 9H X VELe =912F7e4) 
PRINT1215 (V(LoJ)»9 J=1912) 
PRINT1219 (Vtl>sJ)s J=13923) 
121 FORMAT (7/7 9H Y VEL =s1l2F 7-4) 
IF({L~1) 700, 700» 110 
—_ 110 CONTINUE 
CALL DRAW (219X»XNoMOD » O»LAB»s ITITLEse04,5 ONE 0»090969 6sl1sLAST) 
DO 720 N= 2elil 
M=2*N=-1] 
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9 | 
T(2sN)=T( 215M) 

PRINT 77ls Ms U(l29N) sTl29N) 


Ne 


a] 


U(N»,10)= o5/120—.0003 88*UU 
U(Nsl1l)= «3042—-.000230*UU 


GO TO 50 

UU=XM=-149. 

U(Ns2) =1¢4758+-e000290*UU 
U(Ns3) =1-4629+.000266*UU 
U(Ns4) =164239+.000194*UU 
U(N»s5) =123573+.2000092*UU 
U(Ns6) =12¢2611—--000014*UU 
U(N»s7) =121336—-e000104#UU 
U(Ns8) = 09733—e000160*UU 
U(Ns9) = «7796-20001 76*UU 
UCN» 10)= 2«5526—e000150*UU ~ 
U(Nsl1l)= 2«2926—-.000092*UU 
GO TO 50 

UU=XM=-199. 

U(Ns2) =1-4903+e0000121*UU 
U(Ns3) =1¢4762+.0000110*UU 
U(Ns4) =1.24336+.0000080*UU 
U(Ns5) 3=1¢3619+e0000026*UU 
U(INs6) =1-¢2604—-.0000005*UU 
U(Ns7) 


=121284~.0000042*UU 


, ‘ 
; t 


FORMAT (11052F10.5) 
CONTINGE 
L=2 
DO 10 M=29250 
Xi4=M 
Y=O20 
DELY=0el 
~DELX=0.001 
X(M)=DELX*XM 
T=L+1 
N3N=M 
N=L+] 
IF(M= 99) 72973573 
READ 1 (U(ITs11)sU(1510) sU( 139) sU(1I,8) sU( 157) sU( 136) 
1sU( 155) »U( 154) sU6193) sUCITs2) ) 3 
FORMAT (10F604) 
IF(M-99 ) 50951552 
IF(M-149) 515535954 
IF (M-199) 53,;55,55 
UU=XM— 99-6 
U(Ns2) =124388+-000740*#UU 
U(Ns3) =124292+-2000674*UU 
“UC N94) £143993+2000492 UU 
U(Ns5) =1-23454+.000238*UU 
U(N»6) =1.2628—-.000034*UU 
U(Ns7) =1-1467-—-e-000262*UU 
U(Ns8) = 29938—-.000410*UU 
U(Ns9) = «8022—-.000452#UU 
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UINs8) = 269653--0000066*UU 
U(Ns9) = «6 7708—--0000072*UU 
UCN»s10)= 2£5451~-.-0000063*UU 
UC(Ns1l1)= «2880—.0000475*UU 
50 Ui M=O.60 
UM=0.0 
DO 1ll J=2s1l 
XJ=J 
Y(J)SDELY*XJ=_2 
U(L»12)=0-0 
Vi(bL»s2)=0-0 
* V(bsl)=V( 193) 
Y¥(12)=1.0 
T(L9l12)=T(L911)+0e0l 
PCG Py=1 (ls 3) 
U(L»1il)=U(L»3) 
8 VilseJ)= VilLeJ-1) + DELY/(2-*DELX*¥R) * (U(L=15,J)-U(L+1>5J)) 
A= DELX/ (U(LsJ)*#PR*¥DELY**2 ) 
V(lbs12)=0.0 — 
B= (DELX#R¥V(Ls Jy) / (2e*DELY*U(L»J)) 
TIL+tl»sJI=(A-B)¥T( Lot] FOL em2e*¥A)*¥T (Lo) +( AB) ¥T(L Jel) 
IF(J~-2) 13913512 
13 UTM=UTM+e5*U(L» YIAT(L 9 J) 
UM=UM+e5*U( LJ) 
GO TO 14 | 
12 UTM=UTM+U(L»J)*T(LoJ) ” eee 
UM=UM+U(L oJ) aoe 
14 TM=UTM/UM 
mee CONTINUE 
T1= ¢ Tlbh»s12)—-T(L»s1l1)) / DELY 
T2= T(L»1l12)~-TM 
XN(M)=4-0#T1/T2 
IF(N3N=50 ) 315931530 
30 ITF ((N3N/25)*¥25-N3N) 159315931 
31 PRINT 189X(M) oXN(M)oT2eT19TM 
18 FORMAT (// 4H X =9F6e395X914H NUSSELT NOe =9F8e595X98H TW-IM =, 
' 1F82595Xs8H SLOPE =9F8e595X95H TM =9FBe5) 
PRINT 199 (TllLsJ)s J=ls12) 
19 FORMAT (/9H TEMP. =sl2F 764) 
PRINT 209 (U(L9J)s J=Els12) 
20 FORMAT (7/7 9H X VELe =912F 74) 
PRINT 21s (ViLoJ)s J=l912) 
15 DO 60 J=251l 
T(LsJ)HTIL+1l9J) 
U(L~-lsJ)=U(L oJ) 
60 U(LsJ)=U(L4+19J) 
10 CONTINUE i 
IF (KK=6) 76057615761 


761 MOD=3 


GO TO 762 


760 MOD=2 
C 762 CALL DRAW(497 9X »XNoMOD 90 sLABsI TITLE 9004510050» 0» 0906s 691sLAST ), 
750 CONTINUE 


63 


‘ 
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